LINEAR STATIONARY ITERATIVE METHODS FOR THE 
FORCE-BASED QUASICONTINUUM APPROXIMATION 



M. LUSKIN AND C. ORTNER 

Abstract. Force-based multiphysics coupling methods have become popular since 
they provide a simple and efficient coupling mechanism, avoiding the difficulties in 
formulating and implementing a consistent coupling energy. They are also the only 
known pointwise consistent methods for coupling a general atomistic model to a finite 
element continuum model. However, the development of efficient and reliable iterative 
solution methods for the force-based approximation presents a challenge due to the 
non-symmetric and indefinite structure of the linearized force-based quasicontinuum 
approximation, as well as to its unusual stability properties. In this paper, we present 
rigorous numerical analysis and computational experiments to systematically study 
the stability and convergence rate for a variety of linear stationary iterative methods. 



1. Introduction 

Low energy local minima of crystalline atomistic systems are characterized by highly 
localized defects such as vacancies, interstitials, dislocations, cracks, and grain bound- 
aries separated by large regions where the atoms are slightly deformed from a lattice 
structure. The goal of atomistic-to-continuum coupling methods [T141t[T5 | [T6 | l22 | [26 | [28 | 
32] is to approximate a fully atomistic model by maintaining the accuracy of the atom- 
istic model in small neighbors surrounding the localized defects and using the efficiency 
of continuum coarse-grained models in the vast regions that are only mildly deformed 
from a lattice structure. 

Force-based atomistic-to-continuum methods decompose a computational reference 
lattice into an atomistic region A and a continuum region C, and assign forces to rep- 
resentative atoms according to the region they are located in. In the quasicontinuum 
method, the representative atoms are all atoms in the atomistic region and the nodes 
of a finite element approximation in the continuum region. The force-based approxi- 
mation is thus given by 
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where y denotes the positions of the representative atoms which are indexed by j, J^jiy) 
denotes the atomistic force at representative atom j, and ^(y) denotes a continuum 
force at representative atom j. 

The force-based quasicontinuum method (QCF) uses a Cauchy-Born strain energy 
density for the continuum model to achieve a patch test consistent approximation [6j 
[TH 121]. We recall that a patch test consistent atomistic-to-continuum approximation 
exactly reproduces the zero net forces of uniformly strained lattices [19J[211I2Z] • How- 
ever, the recently discovered unusual stability properties of the linearized force-based 
quasicontinuum (QCF) approximation, especially its indefmiteness, present a chal- 
lenge to the development of efficient and reliable iterative methods [12]. Energy-based 
quasicontinuum approximations have many attractive features such as more reliable 
solution methods, but practical patch test consistent, energy-based quasicontinuum 
approximations have yet to be developed for most problems of physical interest, such 
as three-dimensional problems with many-body interaction potentials [20 ] l2T j l30]. 

Rather than attempt an analysis of linear stationary methods for the full nonlinear 
system, in this paper we restrict our focus to the linearization of a one-dimensional 
model problem about the uniform deformation y F and consider linear stationary meth- 
ods of the form 

P{u (n+l) -u {n) ) = ar {n \ (1) 
where P is a nonsingular preconditioning operator, the damping parameter a > is 
fixed throughout the iteration (that is, stationary), and the residual is defined as 

r W := / - Lfu^. 

We will see below that our analysis of this simple model problem already allows us to 
observe many interesting and crucial features of the various methods. For example, we 
can distinguish which iterative methods converge up to the critical strain (see (jSJ) 
for a discussion of the critical strain), and we obtain first results on their convergence 
rates. 

We begin in Sections Eland [3] by introducing the most important quasicontinuum ap- 
proximations and outlining their stability properties, which are mostly straightforward 
generalizations of results from [5HTT|[T5] . In Section HJ we review the basic properties 
of linear stationary iterative methods. 

In Section we give an analysis of the Richardson Iteration (P = I) and prove 
a contraction rate of order 1 — 0(N~ 2 ) in the l v e norm (discrete Sobolev norms are 
defined in Section [27Tj) . where N is the size of the atomistic system. 

In Section [6J we consider the iterative solution with preconditioner P = L^ 1 , where 
Lp Cl is a standard second order elliptic operator, and show that the preconditioned 
iteration with an appropriately chosen damping parameter a is a contraction up to the 
critical strain only in 14 2 '°° among the common discrete Sobolev spaces. We show, 
however, that a rate of contraction in U 2 '°° independent of N can be achieved with the 
elliptic preconditioner L^ cl and an appropriate choice of the damping parameter a. 

In Section [71 we consider the popular ghost force correction iteration (GFC) which 
is given by the preconditioner P = L^ c , and we show that the GFC iteration ceases to 
be a contraction for any norm at strains less than the critical strain. This result and 
others presented in Section [7J imply that the GFC iteration might not always reliably 
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reproduce the stability of the atomistic system [9]. We did not find that the GFC 
method predicted an instability at a reduced strain in our benchmark tests [IB] (see 
also [2l])- To explain this, we note that our ID analysis in this paper can be considered 
a good model for cleavage fracture, but not for the slip instabilities studied in [TEJEI] . 
We are currently attempting to develop a 2D benchmark test for cleavage fracture to 
study the stability of the GFC method. 



2. The QC Approximations and Their Stability 

We give a review of the prototype QC approximations and their stability properties 
in this section. The reader can find more details in l9lfT0l. 



2.1. Function Spaces and Norms. We consider a one-dimensional atomistic chain 
whose 2iV+i atoms have the reference positions Xj = je for e = 1/N. The displacement 
of the boundary atoms will be constrained, so the space of admissible displacements 
will be given by the displacement space 

U={ue R 2N+1 : u_ N = u N = 0}. 

We will use various norms on the space U which are discrete variants of the usual 
Sobolev norms that arise naturally in the analysis of elliptic PDEs. 
For displacements v alA and 1 < p < oo, we define the £ p norms, 



( £ J2e=-N+i M P ) ^ 1 < P < oo, 
max^=_jv+i,...,iv M, p = oo, 



and we denote by U°' p the space U equipped with the £ p norm. The inner product 
associated with the £1 norm is 



N 



(v,w) :— e veVJi for v, w G U. 



-JV+1 



We will also use \\f\\ep and (/, g) to denote the ^-norm and ^-inner product for arbi- 
trary vectors /, g which need not belong to U. In particular, we further define the U ,p 
norm 

hWu 1 * '■= \W\\ef, 

where (v ')<? = v ' t = e~ x {yi — f^-i), £ = —N + 1, . . . , N, and we let U 1,p denote the space 
U equipped with the U 1 ^ norm. Similarly, we define the space U 2 ' p and its associated 
U 2 ' p norm, based on the centered second difference v" = e~ 2 {y£ + i — 2v£ + for 
£ = -N + 1,...,N - 1. 

We have that v' G M 27V for v G U has mean zero Y2j=-N+i v j = ®- ^ e can * nus 
obtain from fL0\ Equation 9] that 

max (u , v') < max (u', a) = \\u\\ui,p < 2 max (u', v'). (2) 
veu ctfr 2JV v&J 
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We denote the space of linear functionals on IA by U* . For g £ U*, s = 0,1, and 
1 < p < oo, we define the negative norms by 

||0||i/-.p := sup (g,v), 

\\v\\ u a,q=l 

where 1 < q < 00 satisfies - + - = 1. We let U~ s,p denote the dual space U* equipped 
with the U~ s ' p norm. 

For a linear mapping A : U\ — > U2 where Ui are vector spaces equipped with the 
norms || • \\u t , we denote the operator norm of A by 

11 a II \\Av\\u 2 

\\A\\l(Ui, U 2 ) ■= SUp —r . 

veu,v^o \\v\\ui 
If U\ = U2, then we use the more concise notation 

H^llwi := ||^4||l(Wi, u±)- 
If A : U ' 2 — > U 0,2 is invertible, then we can define the condition number by 

cond(A) = ||yl||^o,2 ■ ||A _1 ||^o,2. 
When A is symmetric and positive definite, we have that 

cond(A) = X^N-x/Xf 

where the eigenvalues of A are < Xf < ■ ■ ■ < X 2 A N _ 1 . If a linear mapping A : IA — > U 
is symmetric and positive definite, then we define the A-inner product and ^4-norm by 

(v,w) A := (Av,w), \\v\\ 2 A = (Av,v). 

The operator A : 1A\ — > W 2 is operator stable if the operator norm ||^4 _1 ||l(w 2 , Ui) is 
finite, and a sequence of operators Aj : — > is operator stable if the sequence 
|| (^4j) _1 |U(W2 Uij) is uniformly bounded. A symmetric operator A : U 0,2 —> U ' 2 is 
called stable if it is positive definite, and this implies operator stability. A sequence of 
positive definite, symmetric operators Aj : U 0,2 — > U ' 2 is called stable if their smallest 

A 

eigenvalues A x 3 are uniformly bounded away from zero. 

2.2. The atomistic model. We now consider a one- dimensional atomistic chain whose 
2iV + 3 atoms have the reference positions Xj = je for e = 1/N, and interact only with 
their nearest and next-nearest neighbors. 

We denote the deformed positions by i/j, j = —N — 1, . . . , N + 1; and we constrain 
the boundary atoms and their next-nearest neighbors to match the uniformly deformed 
state, = Fje, where F > is a macroscopic strain, that is, 

= ~F(N + l)e, y- N = —FNe, 

Un — FNe, y N+1 = F(N + l)e. 

We introduced the two additional atoms with indices ±(7V + 1) so that y = y F is an 
equilibrium of the atomistic model. The total energy of a deformation y £ IR 27V+3 is 
now given by 

N 



-N 
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where 



N+l 
j=—N 



Vi - Vi-\ 



N+l N+l 
j=—N j=-N+l 



(4) 



Here, is a scaled two-body interatomic potential (for example, the normalized Lennard 
Jones potential, 4>(r) = r~ 12 — 2r~ 6 ), and fj, j = —N, . . . , N, are external forces. The 
equilibrium equations are given by the force balance conditions at the unconstrained 
atoms, 

-J r Hy a ) = f J for j = -N+l,...,N-l, 



y) = Fje for j = —N — 1, -N, N, N + 1, 
where the atomistic force (per lattice spacing s) is given by 
ld£ a (y) 



(5) 



(6) 



We linearize dSJ) by letting w £ IR 27V+3 , w ±A r = u±(n+i) = 0, be a "small" displacement 
from the uniformly deformed state yj = Fje; that is, we define 

= foij = -N-l,...,N+ 1. 

We then linearize the atomistic equilibrium equations (j5]) about the uniformly deformed 
state y F and obtain a linear system for the displacement u a , 

{L F u% = fj for j = -N + l,...,N-l, 

u) = for j = -N- 1, —N, N, N+l, 



where (L F v)j is given by 



(L%v) 



-Vj +1 + 2vj - Vj_i 



Vj +2 + 2Vj - V 



Here and throughout we define 

^:=4>"(F) and ^ := 0"(2F), 

where is the interatomic potential in (j3J). We will always assume that <fi" F > and 
02 F < 0, which holds for typical pair potentials such as the Lennard- Jones potential 
under physically realistic deformations. 

The stability properties of L F can be understood by using a representation derived 
in 0, 

N N 



/ |2 



{L F u, u) = eA f £ 

e=-N+i t=-N 

where A F is the continuum elastic modulus 



e 3 ^ F \u';\ 2 = A F \\u'\\%-e 2 ^ F \W'\ 



2 



(7) 



A F = (f)' F + 4<f>2 F . 
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We can obtain the following result from the argument in [HI Prop. 1] and [12] 
Proposition 1. If 4>2F — 0> then 

(L%U,U) 2 ill 

mm n TTpj — = A F — e v e (p 2F , 

U±N=U±(N+1)=0 



where 



\n"\\ 2 



mm 



U±AT=M±(JV + 1)=0 

satisfies < v e < C for some universal constant C. 

2.2.1. The critical strain F*. The previous result shows that L F is positive definite, 
uniformly as N — > oo, if and only if Ap > 0. For realistic interaction potentials, L F 
is positive definite in a ground state F > 0. For simplicity, we assume that F — 1, 
and we ask how far the system can be "stretched" by applying increasing macroscopic 
strains F until it loses its stability. In the limit as iV — > oo, this happens at the critical 
strain F*, which is the smallest number larger than F , solving the equation 

A Fm = <f>"(F*) + 40" (2F„) = 0. (8) 

2.3. The local QC approximation (QCL). The local quasicontinuum (QCL) ap- 
proximation uses the Cauchy-Born approximation to approximate the nonlocal atom- 
istic model by a local continuum model [5l[23l|26]. For next-nearest neighbor interac- 
tions, the Cauchy-Born approximation reads 

<P (e-\y i+ i - y t -i)) « \ [(j>(2y' e ) + 0(2^ +1 )], 

and results in the QCL energy, for y G IR 2Ar+3 satisfying the boundary conditions ([3]), 

N 

j=-N+l (9) 

<t>{y'-N) + ^v'-n) + <t>(y' N+ i) + ^{iv'n+i) 



+ e 



Imposing the artificial boundary conditions of zero displacement from the uniformly 
deformed state, y^ = Fje, we obtain the QCL equilibrium equations 

-7fy**) = f 3 for J = -N+1,...,N-1, 
yf = Fje for j = -N,N, 



where 



rqcl . A _ ld£^\y) 



n c \y) ■ 



e dy j 



(10) 
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We see from ffTUl) that the QCL equilibrium equations are well-defined with only a single 
constraint at each boundary, and we can restrict our consideration to 
y € R 2Ar+1 with ?/_tv = —F and y^ = F as the boundary conditions. 

Linearizing the QCL equilibrium equations ( jTOl) about y F results in the system 

{Lfu* cl ) 3 = f 3 for j = -N+l,...,N-l, 

uf = for j = —N, TV, 



Lt = A F L 



where 

\f = Ap 

and L is the discrete Laplacian, for v £ U, given by 

" -v j+1 + 2vj - Vj-i 



{Lv)j := -v" 



j 



S 2 



./ v • 1 V 1. (11) 



The QCL operator is a scaled discrete Laplace operator, so 

(Lf l u,u) = Af\\u'\\^2 for all u EU. 

In particular, it follows that L)f l is stable if and only if Ap > 0, that is, if and only if 
F < F*, where F* is the critical strain defined in (jSJ). 

2.4. The force-based QC approximation (QCF). The force-based quasicontin- 
uum (QCF) method combines the accuracy of the atomistic model with the efficiency 
of the QCL approximation by decomposing the computational reference lattice into an 
atomistic region A and a continuum region C, and assigns forces to atoms according to 
the region they are located in. The QCF operator is given by 



7 3 iy) -\Ff(y) ifjGC, (12) 

and the QCF equilibrium equations are given by 

-Jf £ (y«*) = f j for j = -N + l,...,N-l, 

y f = Fje for j = —N, N. 

We note that, since atoms near the boundary belong to C, only one boundary condition 
is required at each end. 

For simplicity, we specify the atomistic and continuum regions as follows. We fix 
K e N, 1 < K < N - 2, and define 

A={-K t ...,K} and C = {-N + 1, . . . , N - 1} \ A. 

Linearizing ( |T2l) about y F , we obtain 

{Lfu^)i = fi for i JV + 1 JV-1, 

uf = for j = -N, N, 

where the linearized force-based operator is given explicitly by 

(L F v)j, for j e A. 



(13) 
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The stability analysis of the QCF operator Lf f is less straightforward [TQIITT] ; we will 
therefore treat it separately and postpone it to Section [3J 

2.5. The original energy-based QC approximation (QCE). The original energy- 
based quasicontinuum (QCE) method [26] defines an energy functional by assigning 
atomistic energy contributions in the atomistic region and continuum energy contribu- 
tions in the continuum region. For our model problem, we obtain 

i>27V+l 



£ qcc (y) = * E s t(v) + £ E £ Xv) for v G 



£€C 



where 



£/(y) = \ (0(2^) + + <f>(y' e+1 ) + 0(2^ +1 )) , and 

£e(y) = KM-i + yd + M) + 4(vt + i) + 0(^+i + vWi)) ■ 

The QCE method is patch tests inconsistent [Ttl8|l2"5|l3T] . which can be seen from the 
existence of "ghost forces" at the interface, that is, VS qce (y F ) = g F ^ 0. Hence, the 
linearization of the QCE equilibrium equations about y F takes the form (see [H Section 
2.4] and Section 2.4] for more detail) 



u 



9? 

qce 



fj 





for j 
for j 



-N + l 
-N, N, 



N 



(14) 



where, for < j < N — 1, we have 



+ 







f 2 Vj - 


- v j- 


1 


F 




e 2 






4- 




+ 2 Vj 


~ V 3 


-2 




As 2 






4- 




+ 2vj 


~ V 3 


-2 




As 2 






4- 


-Vj+2 


+ 2v j 


~ V 3 


-2 




As 2 






4- 


- v j+i 


+ 2v 3 


~ V 3 


-1 




e 2 






4- 




+ 2 Vj 


~ V 3 


-1 




e 2 






4- 




+ 2 Vj 


~ V 3 


-1 




e 2 







+ 



+ 



1 v j+2 - 






e 2e 
2 v j+ i - 


V 3 


1 V j+2 - Vj 


e e 
2 Vj - Vj 


-1 


£ 2£ 
1 Vj - Vj- 2 


£ £ 
1 Vj ~ Vj 


-2 


£ 2£ 


£ 2£ 


• 





< j < K - 2, 
j = K-l, 
J = K, 
, J = K + 1, 
j = K + 2, 
K + 3 < j < N - 1, 



and where the vector of "ghost forces," g, is defined by 



The equations for j 



0, < j < K - 2, 

j = K-l, 
j = K, 
j = K+l, 
j = K + 2, 
0, K + 3 < j < N- 1. 

-N + 1, . . . , — 1 follow from symmetry. 



2^4>2F, 
'2i02F> 
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The following result is a new sharp stability estimate for the QCE operator Lp 20 . Its 
somewhat technical proof is given in Appendix 18.11 

Theorem 2. If K > 1, N > K + 2, and <f)% F < 0, then 

inf (L^ e u, u) = A F + X K (f>2 F , 

IK 11/2=1 

where ~ < \k < 1- Asymptotically, as K -> oo ; we /iave 

A K ~ A, + 0(e~ cK ) where A* « 0.6595 and c » 1.5826. 



2.6. The quasi-nonlocal QC approximation (QNL). The QCF method is the 
simplest idea to circumvent the interface inconsistency of the QCE method, but gives 
non-conservative equilibrium equations [5]. An alternative energy-based approach was 
suggested in [13||33]. which is based on a modification of the energy at the interface. 
The quasi-nonlocal approximation (QNL) is given by the energy functional 

N 

£^\y) := e £ cf>(y' e ) + e £ 4>(y' £ + y' e+1 ) + eJ2i Wv'i) + Mi+i)] > 
i=-N+i ieA i&c 

where we set (j)(y'_ N ) = <fi(y' N+ i) = 0. The QNL approximation is patch test consistent; 
that is, y = y F is an equilibrium of the QNL energy functional. 
The linearization of the QNL equilibrium equations about y F is 

(Lfu^^f) for J = -N+1, 

■ qnl -- for j = —N, N, 



N 



where 



(Lfv 



„ -Vj+i + 2vj - v j_i 



-Vj +2 + 2Vj - Vj-2 

As 2 

—V j+2 + 2Vj - Vj- 2 



+ 



,J j+2 



2v 



As 2 

—v j+ i + 2vj 



Vj-i -Vj + 2Vj_i - Uj_ 2 



0<j<K-l, 
J = K, 
j = K + l, 
2<J<N-1. 



(15) 



We can repeat our stability analysis for the periodic QNL operator in [9j Sec. 3.3] 
verbatim to obtain the following result. 



Proposition 3. If K < N — 1, and <p2F < 0, then 

inf (L^ nl M,w) = A F . 
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Remark 1. Since ft^p = (A F — ft F )/4, the linearized operators (ftp) 1 Lp, (ft' F ) 1 L F Z , 
(ft F )~ l Lf, (ftp)- 1 ^, and (ftp^Lf 1 depend only on A F /ft F , N and K. □ 

3. Stability and Spectrum of the QCF operator 

In this section, we give various properties of the linearized QCF operator, most of 
which are variants of our results in [ini[Tl] . We first give a result for the non-coercivity 
of the QCF operator which lies at the heart of many of the difficulties one encounters 
in analyzing the QCF method. 

Theorem 4 (Theorem 1, If ftp > and ft 2F 6l \ {0} then, for sufficiently 

large N, the operator L^ cf is not positive-definite. More precisely, there exist N G N 
and C 1 >C 2 >0 such that, for all N > N and 2 < K < N/2, 

-dN 1 ' 2 < inf (Lfv,v) < -C 2 N 1 ' 2 . 

IMI,2=1 

The proof of Theorem H] yields also the following asymptotic result on the operator 
norm of Ifp . Its proof is a straightforward extension of [TT1 Lemma 2], which covers 
the case p = 2, and we therefore omit it. 

Lemma 5. Let ft 2 ' F ^ 0, then there exists a constant C3 > such that for sufficiently 
large N, and for 2 < K < N/2, 

r<- 1 M 1 / p < ll/~ qcf ll < C i .,N l / p 

iV S \\L, F \\ L(ul , Py w _i, p) S U 3 iV . 

As a consequence of Theorem |4] and Lemma El we analyzed the stability of L^ cf in 
alternative norms. By following the proof of [101 Theorem 3] verbatim (see also [TO] 
Remark 3]), we can obtain the following sharp stability result. 

Proposition 6. If A F > and ft^p < 0, then L^ cf is invertible with 

\\( L f) 1 \\ L (u°>°°, W^) - 1 / A f- 

If Ap = 0, then L^ cf is singular. 

This result shows that Lp f is operator stable up to the critical strain at which 
the atomistic model loses its stability as well (cf . Section 12 . 2 j) . 

3.1. Spectral properties of L^ cf in U ' 2 = £ 2 . The spectral properties of the L^ cf 
operator are fundamental for the analysis of the performance of iterative methods 
in Hilbert spaces. The basis of our analysis of L^ f in the Hilbert space U 0,2 is the 
surprising observation that, even though L^ f is non- normal, it is nevertheless diago- 
nalizable and its spectrum is identical to that of L^ nl . We first observed this numerically 
in jTOj Section 4.4] for the case of periodic boundary conditions. A proof has since been 
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given in [X3J Section 3], which translates verbatim to the case of Dirichlet boundary 
conditions and yields the following result. 

Lemma 7. For all N > 4, 1 < K < N — 2, we have the identity 

Lf = L~ x Vfh. (16) 

In 'particular, the operator L^ cf is diagonalizable and its spectrum is identical to the 
spectrum of Lf. 

We denote the eigenvalues of L^ nl (and Lf f ) by 

< Xf < ...Xf < ... < A^_i- 

The following lemma gives a lower bound for Xf 11 , an upper bound for A^-i; and 
consequently an upper bound for cond(L^ nl ) = X^^/Xf 11 . 

Lemma 8. If K < N - 1 and <// 2 ' F < , then 

Xf > 2 A F , A^_i < (A F - 4$ F ) e~ 2 = <P" F e-\ and 

cond(Lf ) = ^1 < ( 



A qnl - \y2A 

For the analysis of iterative methods, we are also interested in the condition number 
of a basis of eigenvectors of L^ cf as N tends to infinity. Employing Lemma [7J we 
can write L^ cf = L _1 A qcf L where L is the discrete Laplacian operator and A qcf is 
diagonal. The columns of L~ l are poorly scaled; however, a simple rescaling was 
found in (THl Thm. 3.3] for periodic boundary conditions. The construction and proof 
translate again verbatim to the case of Dirichlet boundary conditions and yield the 
following result (note, in particular, that the main technical step, p31 Lemma 4.6] can 
be applied directly). 

Lemma 9. Let A F > 0, then there exists a matrix V of eigenvectors for the force-based 
QC operator such that cond(V^) is bounded above by a constant that is independent 
ofN. 

3.2. Spectral properties of in U 1,2 . In our analysis below, particularly in Sec- 
tions 16.11 and 16.21 we will see that the preconditioner L^ cl = A F L is a promising 
candidate for the efficient solution of the QCF system. The operator L 1 ! 2 can be un- 
derstood as a basis transformation to an orthonormal basis in U 1,2 . Hence, it will be 
useful to study the spectral properties of L^ cf in that space. The relevant (generalized) 
eigenvalue problem is 

Lfv = XLv, veU, (17) 
which can, equivalently, be written as 

IT 1 Iff v = Xv, veU, (18) 



or as 



L~ l ' 2 LfL^ 2 w = Xw, we U, (19) 
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with the basis transform w = L l / 2 v, in either case reducing it to a standard eigenvalue 
problem in £ 2 . Since L and L 1 ! 2 commute, Lemma [7] immediately yields the following 
result. 

Lemma 10. For all N > 4, 1 < K < N — 2 the operator L~ 1 L ( p i is diagonalizable 
and its spectrum is identical to the spectrum of L^Lp^ 1 . 

We gave a proof in [12] of the following lemma, which completely characterizes 
the spectrum of L" 1 /^ 111 , and thereby also the spectrum of L~ 1 L F zf . We denote the 
spectrum of L~ l Lf l (and L~ l Lf) by {/if 1 : j = 1, . . . , 2N - 1}. ' 

Lemma 11. Let K < N — 2 and Ap > 0, then the (unordered) spectrum of L" 1 ^* 11 
(that is, the IA 1 ' 2 -spectrum) is given by 

qnl= fA F -40 2 Vsin 2 ( i ^ i ), j = l,...,2K+l, 
^ ~\A F , j = 2K + 2,...,2N-l. 

In particular, if 0' 2 ' F < 0, then 

max i ^ qnl , 40' 2 ' F . 2 f i 2K+i]n\ 4>" F , M'f • 2/ tt \ ,nrr-2\ 

— — ij = 1 -^— sm = r + i~ sm u^Ti) = + )• 

millj fj,j Ap \ / Ap Ap Ap 

We conclude this study by stating a result on the condition number of the matrix 
of eigenvectors for the eigenvalue problem ffl9l) . Letting V be an orthogonal matrix of 
eigenvectors of L~ l ^ 2 L F al L~ 1 ^ 2 and A the corresponding diagonal matrix, then Lemma 
[7J yields 

L-^LfL- 1 ' 2 = L-^L-^LfL-^L 
= (V T L)- l ~A(V T L). 

Clearly, cond(V" T L) = 0(N 2 ), which gives the following result. 

Lemma 12. If Ap > 0, then there exists a matrix W of eigenvectors for the 
preconditioned force-based QC operator L~ 1 / 2 L F Z L~ x l 2 , such that cond(W) = 0(N 2 ) 
as N oo. 

4. Linear Stationary Iterative Methods 

In this section, we investigate linear stationary iterative methods to solve the lin- 
earized QCF equations ( TT5j) . These are iterations of the form 

P{u {n) -u {n ~ l) ) =ar {n ~ l \ (20) 

where P is a nonsingular preconditioner, the step size parameter a > is constant 
(that is, stationary), and the residual is defined as 

r (n) := / - Lfu {n) . 

The iteration error 

e (n) . = ^qcf _ u (n) 
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satisfies the recursion 
or equivalently, 



Pe (n) = (p - aLfy^, 



e W = (/ - aP- l Lf)e^ =: Ge^, (21) 

where the operator G = I — aP~ l L ( f: i : U — > U is called the iteration matrix. By 
iterating (12"T|) . we obtain that 

e (») = (/ - aP^LfYeW = G n e®. (22) 

Before we investigate various preconditioners, we briefly review the classical theory 
of linear stationary iterative methods [29]. We see from (|22|) that the iterative method 
( 12"U|) converges for every initial guess E IA if and only if G™ — > as n — > oo. For a 
given norm \\v\\, for v G W, we can see from f[2"2l that the reduction in the error after 
n iterations is bounded above by 

lle^ll 

||G n || = sup || . . . 

It can be shown [29] that the convergence of the iteration for every initial guess 
G W is equivalent to the condition p(G) < 1, where p(G) is the spectral radius of 

G, 

p(G) = max { | Aj| : Aj is an eigenvalue of G} . 
In fact, the Spectral Radius Theorem [29] states that 

lim \\G n \\ 1/n = p{G) 



for any vector norm on U. However, if p(G) < 1 and ||G|| > 1, the Spectral Radius 
Theorem does not give any information about how large n must be to obtain < 1. 

On the other hand, if p{G) < 1, then there exists a norm || • || such that ||G|| < 1, 
so that G itself is a contraction [17]. In this case, we have the stronger contraction 
property that 

||e (n) || < llGlllle^H < ||G|n|e (0) ||. 
In the remainder of this section, we will analyze the norm of the iteration matrix, 
||G||, for several preconditioners P, using appropriate norms in each case. 

5. The Richardson Iteration (P = I) 

The simplest example of a linear iterative method is the Richardson iteration, where 
P — I. If follows from Lemma M that there exists a similarity transform S such that 

Lf = S^A^S, (23) 

where cond(S') < C (where C is independent of N), and A qnl is the diagonal matrix 
of W°' 2 -eigenvalues (A| nl )|^f 1 of Lf { . As an immediate consequence, we obtain the 
identity 

G id (a) = I - aLf = S' 1 (i - aA qnl )S, 

where yields 

||Gid(a)||^ < cond(S)||/ - aA qnl || £i < C _max ^ |l - aXf\. (24) 
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If A F > 0, then it follows from Proposition [3] that AJ 11 > for all j, and hence that 
the iteration matrix Gid(a) .— I — aL^ f is a contraction in the || ■ ||^ norm if and only 
if < a < a\f mx := 2/Ag_ 1 . It follows from Lemma EJthat < ax < (2e 2 )/<p F . 

We can minimize the contraction constant for G;d(a) in the \\v\\ $ts norm by choosing 
a = ajjpt := 2/(A ( ] inl + X^n-i)^ an< ^ m this case we obtain from Lemma |8]that 



It thus follows that the contraction constant for Gid(a) in the norm is only of the 

order 1 — 0(e 2 ), even with an optimal choice of a. This is the same generic behavior 
that is typically observed for Richardson iterations for discretized second-order elliptic 
differential operators. 

5.1. Numerical example for the Richardson Iteration. In Figure [H we plot the 
error in the Richardson iteration against the iteration number. As a typical example, 
we use the right-hand side 

{1, x > 0, 
_ ' (25) 
—1, x < 0, 

which is smooth in the continuum region but has a discontinuity in the atomistic 
region. We choose (p" F = 1, Ap = 0.5, and the optimal a = a^ t discussed above (we 
note that Gid(o4p t ) depends only on Af/4>" f and N, but depends on Ap and 4> F 
independently) . We observe initially a much faster convergence rate than the one 
predicted because the initial residual for ( 125]) has a large component in the eigenspaces 
corresponding to the intermediate eigenvalues A^ nl for 1 < j < 2N — 1. However, after 
a few iterations the convergence behavior approximates the predicted rate. 



6. Preconditioning with QCL (P = L^ cl = A F L) 

We have seen in Section [5] that the Richardson iteration with the trivial precondi- 
tioner P = I converges slowly, and with a contraction rate of the order 1 — 0(e 2 ). The 
goal of a (quasi-) optimal preconditioner for large systems is to obtain a performance 
that is independent of the system size. We will show in the present section that the 
preconditioner P = ApL (the system matrix for the QCL method) has this desirable 
quality. 

Of course, preconditioning with P = ApL comes at the cost of solving a large linear 
system at each iteration. However, the QCL operator is a standard elliptic operator 
for which efficient solution methods exist. For example, the preconditioner P = ApL 
could be replaced by a small number of multigrid iterations, which would lead to a 
solver with optimal complexity. Here, we will ignore these additional complications 
and assume that P is inverted exactly. 

Throughout the present section, the iteration matrix is given by 

G qd (a) := / - aiLfY^f = 1- a(A F L)~ l Lf { , (26) 
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Figure 1. Normalized £j?-error of successive Richardson iterations for 
the linear QCF system with N = 200, K = 8, 32, (j)" F = 1, A F = 0.5, 
optimal a = a^ pt , right-hand side fl25]) . and starting guess = 0. 

where a > and Ap — 4>"f + ^4>2F > 0- We will investigate whether, if IA is equipped 
with a suitable topology, G qc \(a) becomes a contraction. To demonstrate that this is 
a non-trivial question, we first show that in the spaces 1 < p < oo, which are 
natural choices for elliptic operators, this result does not hold. 

Proposition 13. If 2 < K < N/2, (p'^p ^ 0, and p G [1, oo), then for any a > we 

have 

||Gqci(a)|| wllP -N 1 ^ asN^oo. 
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Proof. We have from ([2]) and q = p/(p — 1) the inequality 
||L _1 L^ cf || wlp = max || (L^L^u)' 11 



\\ u '\\eP =1 

<2 max ( (L^L^uY, v 

||u'|Lp=l, ||i/|L?=l 

= 2 max (l{L~ 1 L^u),v 

||u'|Lp=l, ||f/|Lg=l 

= 2 max (ip f u, 

||u'|| £ p=l, ||u'|| £ 8=l 



qcf I 



as well as the reverse inequality 



-qcf II ^ II r -1 rqcf II 
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The result now follows from the definition of G qc \(a) in (126]) . Lemma [51 and the fact 
that a > and A F > 0. □ 

We will return to an analysis of the QCL preconditioner in the space U 1,2 in Section 



16. 3[ but will first attempt to prove convergence results in alternative norms. 

6.1. Analysis of the QCL preconditioner in U 2,oc . We have found in our previous 
analyses of the QCF method [iniELTI that it has superior properties in the function 
spaces U 1,0 ° and U 2 '°°. Hence, we will now investigate whether a can be chosen such 
that Gqci(a) is a contraction, uniformly as N — > oo. In [10], we have found that 
the analysis is easiest with the somewhat unusual choice IA 2 '°°. Hence we begin by 
analyzing G qc \(a) in this space. 

To begin, we formulate a lemma in which we compute the operator norm of G qc i(a) 
explicitly. Its proof is slightly technical and is therefore postponed to Appendix 18.21 



Lemma 14. If N > A, then 

ll^qcl(«)|| W 2,o 



1 -<*(!-*) 



a 



2 <t>2F 



A F 



What is remarkable (though not necessarily surprising) about this result is that the 
operator norm of G qc \(a) is independent of iV and K. This immediately puts us into 
a position where we can obtain contraction properties of the iteration matrix G qc i(c>!), 
that are uniform in N and K. It is worth noting, though, that the optimal contraction 
rate is not uniform as Af approaches zero; that is, the preconditioner does not give 
uniform efficiency as the system approaches its stability limit. 

Theorem 15. Suppose that N > 4, A F > 0, and 0' 2 ' F < 0, and define 

qcl,2,oo _ A_F _ 2A F , qcl,2,oo _ ^Ap 

a ° pt - a f + 2\<% f \- y F + A F " max 

Then G qc i(aO is a contraction of U 2,oc if and only if < a < a^ 2 ' 00 , and for any 
such choice the contraction rate is independent of N and K . The optimal choice is 

,qcl,i 
opt 



a = <y2pi 2, °° > which gives the contraction rate 



1_ A F 

( qcl,2,oo\|| _ ~¥f 

I'jqcl^opt )\\U 2 >°° ~ 1+ A F ^ l - 



Proof. Note that «opt 2 '°° = — Hence, if we assume, first, that < a < 



qcl.2,oo ,i 

"opt > then 



\\G qcl (a)\\ u ^ = l- a (l- 2§f) - 2a§f = 1 - a =: m x {a). 
The optimal choice is clearly a = ct^ 2 ' 00 which gives the contraction rate 



I / qcl,2,oo\ II _ qcl,2,oo 

l^qcl^opt J||^2,oo "opt 



20'' 



2F 



A, 
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Alternatively, if a > a^f ,0 °, then 

llGqdCaOHaa.oo = «(l - ^p) - 1 = «^ - 1 = : m 2 (a). 

This value is strictly increasing with a, hence the optimal choice is again a = a2pt 2 '°°- 
Moreover, we have 7712(0;) < 1 if and only if 

a< ^L = a qci,2,oc 

ill "max 

Since, for a = a^ 2 ' 00 we have Tn\(a) = m 2 (a) < 1, it follows that a-mlx' 00 > a opt 2,0 ° 
(as a matter of fact, the condition a^lx' 00 > ^opt' 2 ' 00 is equivalent to Ap > 0). In 
conclusion, we have shown that ||G qc i(a)||^2,oo is independent of N and K and that it 
is strictly less than one if and only if a < a^^' 00 , with optimal value a = aT^ 2 ' 00 . □ 

As an immediate corollary, we obtain the following general convergence result. 

Corollary 16. Suppose that N > 4, Ap > 0, ifi'^p < 0, and suppose that || • \\x is a 
norm defined on U such that 

\\ u \\x < C|| M llw 2 .°° Vu G IA. 

Moreover, suppose that < a < a^'°°. Then, for any u ElA, 

\\G qc i{ct) n u\\ x < (7 n C||w||^2,oo — >• as n — >• 00, 

where q := \\G qc i(a)\\u^.°° < 1. 

In particular, the convergence is uniform among all N, K and all possible initial 
values u ElA for which a uniform bound holds. 

Proof. We simply note that, according to Theorem [T5l for < a < a^x'°°; we have 

||Gqd(a) n || Ma .cc < <f\ 

where q := ||C qc i(ct) ||^a,oo < 1 is a number that is independent of N and K. Hence, we 
have 

||G qc i(oO n w|| x < C\\G qcl {a) n u\\ u2>oa < Cq n \\u\\ U 2, x . 

□ 



Remark 2. Although we have seen in Theorem US] and Corollary [TO] that the linear 
stationary method with preconditioner ApL and with sufficiently small step size a 
is convergent, this convergence may still be quite slow if the initial data is "rough." 
Particularly in the context of defects, we may, for example, be interested in the con- 
vergence properties of this iteration when the initial residual is small or moderate in 
IA 1,P , for some p E [1, 00], but possibly of order O(N) in the W 2, °°-norm. We can see 
from the following Poincare and inverse inequalities 

IMIw 1 ' 00 < ^IMIw 2 - 00 an d |l M llw 2 .°° < 27V||ii|| w i >0 o for all u ElA] 
that the application of Corollary [16] to the case X = U 1 ' 00 gives the estimate 
||G f qoi(a!) n M|| W i,oo < q n N\\u\\ u i,°o for all u E U. 
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Similarly, with X = U 1,2 , we obtain 

\\G qcl (a) n u\\ ul>2 < q n N 3/2 \\u\\ u i,2 for all u^U. (27) 

We have seen in Proposition [TBI that a direct convergence analysis in W 1,p , p < oo, 
may be difficult with analytical methods, hence we focus in the next section on the 
case U 1 ' 00 . □ 

6.2. Analysis of the QCL preconditioner in lA 1,oa . As before, we first compute 
the operator norm of the iteration matrix explicitly. The proof of the following lemma 
is again postponed to the Appendix 18.21 



Lemma 17. // K > 3, N > max(9, K + 3), and <p'^ F < 0, then 

fa 

<1-2S)I + «(6 + ^ ~ /or ag 1 ' 00 < a 



|G qc i(a)|| W i 



a\+aA\^f\ /or < a < og 1 ' 00 , 



l + (2 + e-2eA)fe| 



qcl.l.oo 
rv ■ = 

"opt 

i- r> qcl,2,oo ^ qcl.l.oo ^ -i 

satisfies a£ pt ' < a£ pt ' < 1 . 

Again we note that the operator norm is independent, but now up to terms of order 
0(eK), of the system size. 

Theorem 18. Suppose that K > 3, N > max(9, K + 3) , and 0^ < 0, i/ien £/ie 
following statements are true: 

(i) If (f>p + 8(f>2F < 0, £/ien G qc i(a) is not a contraction ofU 1,oc , for any value of a. 

(ii) If 4> F + 8<p2F > 0, then G qc i(a) is a contraction for sufficiently small a. More 
precisely, setting 

max • A F + (8 + 2e-4eK)\<% F \' 

we have that G qc i(o;) is a contraction of U l,oc if and only if < a < o^ax' 00 - 
The operator norm \\G qc i(a)\\i(i,oo is minimized by choosing a = a™^ 1 ' 00 (cf. 
Lemma\Tl^ and in this case 

M^y / qcl,l,oo\|| _ 'Pf + , -I 

l^qcl^opt ) || M l,oc - ^ + (2 _ £ + 2ejK -)^ < ^ 

Proof. Suppose, first, that < a < a^t' 1 ' 00 . Since ag 1 ' 00 < 1 it follows that 

|p qc i(a)|| liDO = 1 - " 2 ' 

Ap 

and hence ||G qc i(a)||^i,°o < 1 if and only if <p" F + 80 2 ' F > 0- m that case ||G qc i(a)||^i,oo 
is strictly decreasing in (0, ag 1 ' 00 ]. 

Since ag 1 ' 00 > a^f' 00 = (1— 2^) _1 we can see that ||G qc i(a)|| w i,oo is always strictly 
increasing in [ojopt' 1 ' 00 , +oo) and hence if <f) F + 80 2 V > 0, then a = o^ 1,oc minimizes 
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the operator norm ||G qc i(a)||^i,°o. Moreover, straightforward computations show that 
"mi 1 ' 00 > "opt 1 ' 00 and that ||G qd (a)|| M i,oo < 1 if and only if < a < a^'°°. □ 

We remark that the optimal value of a in W 1,0 °, that is a = ttopt' 1 ' 00 ; is n °t the 
same as the optimal value, a^ 2 ' 00 , i n U 2, °° . However, it is easy to see that a^ l, °° = 
"opt' 2 ' 00 + 0(eK), and hence, even though a^t' 2 ' 00 is not optimal in U l,oa it is still close 
to the optimal value. On the other hand, a^ 1,ra and a*^ 2 ' 00 are not close, since, if 
AeK - 2e < 1, then 

qcl,l,oo < 2A F ZA F _ qc l >2 ,oo 

max -^ + 3|0^l 0f ^ ' 

In summary, we have seen that the contraction property of G qc i(a) in U l,oc is signif- 
icantly more complicated than in W 2 '°°, and that, in fact, G qc \(a) is not a contraction 
for all macroscopic strains F up to the critical strain F*. 



6.3. Analysis of the QCL preconditioner in U 1,2 . Even though we were able 
to prove uniform contraction properties for the QCL-preconditioned iterative method 
in W 2, °°, we have argued above that these are not entirely satisfactory in the pres- 
ence of irregular solutions containing defects. Hence we analyzed the iteration matrix 
G qc \(a) = I — a(ApL)~ 1 L c ^ i in U 1,0 °, but there we showed that it is not a contraction 
up to the critical load F*. To conclude our results for the QCL preconditioner, we 
present a discussion of G qc \(a) in the space U 1,2 . 
We begin by noting that it follows from ( 1211) that 

pl/2 e (n) = pl/2 Gqd(a)e (n-l) = pl/2 (7 _ a p-l L f\ p-1/2 (pV2 e (n-l)) 

= (/ - aP-^Lf P" 1 / 2 ) (P^e^) =: G qcl (a) (P^V*" 1 )) . 

Since HP 1 / 2 !; \\p = A]/ 2 ||f ||^i,2 for v E U, it follows that G qc i(a) is a contraction in 
U 1,2 if and only if G qc \(a) is a contraction in £ 2 . Unfortunately, we have shown in 
Proposition [TBI that ||G qc i(a)||^i,2 ~ N 1 / 2 as N — > oo. Hence, we will follow the idea 

used in Section and try to find an alternative norm with respect to which G qc \(a) is 
a contraction. 

From Lemma [TU] we deduce that there exists a similarity transform S such that 
cond(iS) < iV 2 , and such that 

L~ l ' 2 LfL~ 1 ' 2 = S^A^S, 

where A qnl is the diagonal matrix of ZY^-eigenvalues (/x^ 111 ) 2 ^ -1 of L]? 1 . As an imme- 
diate consequence we obtain 



G qcl (a) = S" 1 (l-^A^)S. 
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Proceeding as in Section [5J we would obtain that ||G qc i(a)||£2 < 0(N 2 ). Instead, we 
observe that 



G , qc i(a)ML T c = ||SGqci(a)w|| 2 = 11(7- f-A qnl )Su 



that is, 



\G^(a)\\ §T§ < , = max _ i |1 - ^nf\. (28) 



Thus, we can conclude that G qc \(a) is a contraction in the || ■ II^T^-norm if and only if 
< a < a^zx' 2 := / fJ>2N-v Moreover, we obtain the error bound 

||e (n) || w i, 2 < cond(S) g n ||e (0) || w i, 2 < N 2 q n \\e {0) || w i,a, 



where q := ||G qc i(a) \\g T g- This is slightly worse in fact, than ( )27|) . however, we note 
that this large prefactor cannot be seen in the following numerical experiment. 

Moreover, optimizing the contraction rate with respect to a leads to the choice 
a opt' 1,2 := ^Ap/i^" 1 + A^tv-i)' an d m this case we obtain from Lemma [TT] that 

qnl qnl \ _ Af 

r, - „ —\\r (n,^ d ^ 2 \\\ - /i2Af ~ 1 ZJ^I < #L 

H Vopt • ||^qclV«opt )\\sTg qnl qnl — 1 , Ap ' 

V2N-1 + A*i ^ 4> F 

where the upper bound is sharp in the limit K — > 00. It is particularly interesting to 
note that the contraction rate obtained here is precisely the same as the one in U 2 '°° (cf. 
Theorem [T5|) . Moreover, it can be easily seen from Lemma [TTl that a^' 1 ' 2 — > ckopt 2 '°° 
as K — > 00, which is the optimal stepsize according to Theorem [T51 We further have 
that a&£> 2 -> a^°° as K ->• 00. 

IIleLX IIleLX 

6.4. Numerical example for QCL-preconditioning. We now apply the QCL- 
preconditioned stationary iterative method to the QCF system with right-hand side 
( )2"5|) . 4>"f = 1' Af = 0.2, and the optimal value a = c^' 2 ' 00 (we note that Gid{a^ 2 '°°) 
depends only on Ap/(j)" F and N, but depends on A F and (p" F independently). The 
error for successive iterations in the U 1,2 , U 1,oc and W 2, °°-norms are displayed in Figure 
HJ Even though our theory, in this case, predicts a perfect contractive behavior only 
in U 2,co and (partially) in U 1 ' 2 , we nevertheless observe perfect agreement with the 
optimal predicted rate also in the W 1,00 -norms. As a matter of fact, the parameters 
are chosen so that case (i) of Theorem [18] holds, that is, G qc \(a) is not a contraction of 
U l,oc . A possible explanation why we still observe this perfect asymptotic behavior is 
that the norm of G qc \(a) is attained in a subspace that is never entered in this iterative 
process. This is also supported by the fact that the exact solution is uniformly bounded 
in U 2,oc as N, K — > 00, which is a simple consequence of Proposition El 

7. Preconditioning with QCE (P = L^ c ): Ghost-Force Correction 

We have shown in [5l[T2] that the popular ghost force correction method (GFC) is 
equivalent to preconditioning the QCF equilibrium equations by the QCE equilibrium 
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Figure 2. Error of the QCL-preconditioned linear stationary iterative 
method for the QCF system with N = 800, K = 32, <p F = 1, A F = 0.2, 
optimal value a = «opt 2 '°°, and right-hand side (1231) . In this case, the 
iteration matrix G qc \(a) is not a contraction of U 1,oc . Even though our 
theory predicts a perfect contractive behavior only in U 2 '°°, we observe 
perfect agreement with the optimal predicted rate also in the U 1 ' 2 and 
W 1,oc -norms. 

equations. The ghost force correction method in a quasi-static loading can thus be 
reduced to the question whether the iteration matrix 

r< T ( r QCCN — 1 rqcf 

<J"qce •— 1 — V-I^F ) F 

is a contraction. Due to the typical usage of the preconditioner L^ CG in this case, we 
do not consider a step size a in this section. The purpose of the present section is (i) 
to investigate whether there exist function spaces in which G qcc is a contraction; and 
(ii) to identify the range of the macroscopic strains F where G qcc is a contraction. 
We begin by recalling the fundamental stability result for the L^ ce operator, Theorem 

m 

inf (Lf c u, u) = A F + \k<I>2Fi 

ll«'ll £ |=l 

where Xk ~ A* + 0(e~ cK ) with A* ~ 0.6595. This result shows that the GFC iteration 
must necessarily run into instabilities before the deformation reaches the critical strain 
F*. This is made precise in the following corollary which states that there is no norm 
with respect to which G qcc is a contraction up to the critical strain F*. 

Corollary 19. Fix N and K , and let || • \\x be an arbitrary norm on the space U, 
then, upon understanding G qcc as dependent on <j)" F and 4>2F> we have 
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Despite this negative result, we may still be interested in the question of whether the 
GFC iteration is a contraction in "very stable regimes," that is, for macroscopic strains 
which are far away from the critical strain F*. Naturally, we are particularly interested 
in the behavior as N — > oo, that is, we will investigate in which function spaces the 
operator norm of G qce remains bounded away from one as iV — > oo. Theorem H] on the 
unboundedness of L^ cf immediately provides us with the following negative answer. 

Proposition 20. If 2 < K < N/2, ^ F ^ 0, and A F + X K (fi% F > 0, then 

\\G qce \\ui,2 ~ N 1/2 as N -> oo. 

Proof. It is an easy exercise to show that, if Ap + Xk^f > 0' then the W 1,2 -norm is 
equivalent to the norm induced by L^ e , that is, 

C _1 ||m|| w i,2 < \\u\\ L v» < C||W|| W 1,2. 

Hence, we have ||Cr qce ||^i,2 ~ IjGqcelljq^ and by the same argument as in the proof of 
Proposition [T3J and using again the uniform norm-equivalence, we can deduce that 

|| G qce|| w i,a ~ H^L^La, U~^) ± 1 ~ N ^ as N y OO. 

□ 

Since the operator (L^ 00 )" 1 /^ is more complicated than that of (ApL)~ 1 L F 2t , which 
we analyzed in the previous section, we continue to investigate the contraction proper- 
ties of G qce in various different norms in numerical experiments. In Figure [3], we plot 
the operator norm of G qce , in the function spaces 

U k ' p , fc = 0,l,2, p= l,2,oo, 

against the system size TV (see Appendix 18.31 for a description of how we compute 
||G qce ||^fc,p). This experiment is performed for Ap/^'p = 0.8 which is at some distance 
from the singularity of L^ ce (we note that G qce depends only on Ap/<p F and N since 
both (<p" F )~ l Lp 21 and (<p" F )~ l Lnp C depend only on Ap/(p F and N). The experiments 
suggests clearly that ||Gq C o||wfe p — > 00 as N — > 00 for all norms except for U 1 ' 00 and 
U 2 ' 1 . 

Hence, in a second experiment, we investigate how ||Grq Ce ||wi,°° and ||G qce ||w2,i be- 
have, for fixed N and K, as Ap + Xk4>2f approaches zero. The results of this exper- 
iment, which are are displayed in Figure HI confirm the prediction of Corollary [19] 
that ||G qcc || W fc,p — > 00 as A F + Xr^f approaches zero. Indeed, they show that 
IIGqcell^fe.p > 1 already much earlier, namely around a strain F where Ap w 0.52 
and Ap + Xk^f ~ °- 44 - 

Our conclusion based on these analytical results and numerical experiments is that 
the GFC method is not universally reliable near the limit strain F*, that is, under 
conditions near the formation or movement of a defect it can fail to converge to a 
stable solution of the QCF equilibrium equations as the quasi-static loading step tends 
to zero or the number of GFC iterations tends to infinity. Even though the simple model 
problem that we investigated here cannot, of course, provide a definite statement, it 
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Operator Norm of G qce = I — (ij ce ) 1 L'^ i 




Figure 3. Graphs of the operator norm ||G qce ||^fe,p, k — 0,1,2, p — 
1,2, oo, plotted against the number of atoms, N, with atomistic region 
size K = \V~N~\ — 1, and Ap/cf)" F = 0.8. (The graph for the W^-norms, 
p = l,oo, are only estimates up to a factor of 1/2; cf. Appendix 18.31 ) 
The graphs clearly indicate that ||G qce ||^fe,p — > oo as N — > oo in all spaces 
except for U 1 ' 00 and U 2 ' 1 . 



shows at the very least that further investigations for more realistic model problems 
are required. 



Conclusion 

We proposed and studied linear stationary iterative solution methods for the QCF 
method with the goal of identifying iterative schemes that are efficient and reliable for 
all applied loads. We showed that, if the local QC operator is taken as the precondi- 
tioned then the iteration is guaranteed to converge to the solution of the QCF system, 
up to the critical strain. What is interesting is that the choice of function space plays a 
crucial role in the efficiency of the iterative method. In U 2 '°°, the convergence is always 
uniform in N and K, however, in U 1,0 ° this is only true if the macroscopic strain is at 
some distance from the critical strain. This indicates that, in the presence of defects 
(that is, non-smooth solutions), the efficiency of a QCL-preconditioned method may 
be reduced. Further investigations for more realistic model problems are required to 
shed light on this issue. 

We also showed that the popular GFC iteration must necessarily run into instabilities 
before the deformation reaches the critical strain F*. Even for macroscopic strains that 
are far lower than the critical strain F*, we show that ||G qce ||^i,2 ~ N 1 ^ 2 . We then give 
numerical experiments that suggest that ||G qcc ||^fe, P — > oo as N — > oo for all tested 
norms except for U 1 ' 00 and U 2 ' 1 . 
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The results presented in this paper demonstrate the challenge for the development of 
reliable and efficient iterative methods for force-based approximation methods. Further 
analysis and numerical experiments for two and three dimensional problems are needed 
to more fully assess the implications of the results in this paper for realistic materials 
applications. 




FIGURE 4. Graphs of the operator norm ||G qce || W fe, P , (k,p) G 
{(l,oo),(2,l)}, for fixed N = 256, K = 15, <\>" F = 1, plotted against 
Ap. For the case U 1,co only estimates are available and upper and lower 
bounds are shown instead (cf. Appendix 18. 3j) . The graphs confirm the 
result of Corollary [191 that ||Cr qC e||wfc,p ->ooas Ap + ^k4>2F ~~ >" 0- More- 
over, they clearly indicate that ||G qce ||^fc,p > 1 already for strains F in 
the region Ap ~ 0.5, which are much lower than the critical strain at 
which L^ ce becomes singular. 
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8. Appendix 

8.1. Proof of Theorem [2l The purpose of this appendix is to prove the sharp stability 
result for the operator Lf°, formulated in Theorem EJ Using Formula (23) in [9] we 
obtain the following representation of L^, 

-K-2 N 



(Lp Ce n,u)=< 2^ eA F \u' e \ 2 + ^ £Af|i4| 



-JV+l l=K+3 
K-l 



£ e(A F \u' e \ 2 -e 2 ^M\ 2 ) 
e=-K+2 ) (29) 

+ e{(A F - ^ F )(K K+1 \ 2 + K| 2 ) + A F (\u'_ K \ 2 + \u' K+1 \ 2 ) 

+ (A F + ^ F )(\u'_ K _ 1 \ 2 + \u' K+2 \ 2 ) 

1 ^2 ill n II \2 i I // 1 2 i I II \2 i I // |2\ 1 

-2 £ ^2f(I«-xI +Kk-iI + \u k \ + \u K+1 \)y 

If 4>2 F < 0, then we can see from this decomposition that there is a loss of stability 
at the interaction between atoms —K — 2 and —K — 1 as well as between atoms K + 1 
and K + 2. It is therefore natural to test this expression with a displacement u defined 
by 



fl'r 




From ( 1291) . we easily obtain 



{Lf%u) = A F + yi 



2F- 



In particular, we see that, if A F + \<p2 F < 0, then Up e is indefinite. On the other 
hand, it was shown in [8] that Lp 20 is positive definite provided A F + 2 'f > 0- a 
matter of fact, the analysis in [8] is for periodic boundary conditions, however, since 
the Dirichlet displacement space is contained in the periodic displacement space the 
result is also valid for the present case.) 
Thus, we have shown that 

inf (Lf e u, u) = A F + jU02F; where | < fj, < 1. 

IKIIf2=l 

To conclude the proof of Theorem [2], we need to show that /i depends only on K and 
that the stated asymptotic result holds. 

From (129]) it follows that Up e can be written in the form 

{L^u^u) = (u') T Uu\ 
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where we identify u' with the vector vl 



u 



e)e=-N+i 



and where H G 



t>2Nx2N 



Writing 



H = 4>"pHi + 4>'npH.ii w e can see that Hi = Id and that "H 2 has the entries 



/ 



U 2 



\ 



l 

2 1 
1 3/2 
1/2 



1/2 
3 1/2 
1/2 9/2 




v ■■/ 

Here, the row with entries [1, 3/2, 1/2] denotes the Kth row (in the coordinates u' k ). 
This form can be verified, for example, by appealing to fT29|) . Let cr(A) denote the 
spectrum of a matrix A. Since, by assumption, 2F < 0, the smallest eigenvalue of Ji 
is given by 

min (t(H) = 4>'p + (\)' 2F max (riHo), 

that is, we need to compute the largest eigenvalue A of "H 2 - Since = for 

k = K + 3, K + 4, . . . and for K = —K — 2, —K — 3, . . . , and since eigenvectors 
are orthogonal, we conclude that all other eigenvectors depend only on the submatrix 
describing the atomistic region and the interface. In particular, A depends only on K 
but not on N. This proves the claim of Theorem [2] that Xk depends indeed only on K. 



We thus consider the {—K — 1, 



/ 9/2 
1/2 



1/2 
3 
1/2 



1/2 
3/2 
1 



K + 2}-submatrix "H 2; which has the form 

\ 



V 



1 

3/2 
1/2 



1/2 
3 
1/2 



1/2 
9/2 / 



-K + 2, 



K 



Xijji 



Letting H 2 ^ = A^, then for £ 

^ t _ x + 2rj) t + ipi+i 
and hence, ip has the general form 

= az e + bz~ l , e=-K+l, 
1,K + 1,K + 



leaving ipe undefined for £ £ {—K,—K 
are the two roots of the polynomial 



...,K, 

2} for now, and where z,l/z 



+ (2-A)z + l = 0. 



In particular, we have 



1 > 1. 



(30) 
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To determine the remaining degrees of freedom, we could now insert this general 
form into the eigenvalue equation and attempt to solve the resulting problem. This 
leads to a complicated system which we will try to simplify. 

We first note that, for any eigenvector ip, the vector (ipK-e) is also an eigenvector, 
and hence we can assume without loss of generality that ip is skew-symmetric about 
£ = 1/2. This implies that a = —b. Since the scaling is irrelevant for the eigenvalue 
problem, we therefore make the ansatz ipi = z e — z~ l . Next, we notice that for K 
sufficiently large the term z~ l is exponentially small and therefore does not contribute 
to the eigenvalue equation near the right interface. We may safely ignore it if we 
are only interested in the asymptotics of the eigenvalue A as K — > oo. Thus, letting 
ipt — z , £ = 1, . . . , K and ipe unknown, £ = K + 1, K + 2, we obtain the system 



- + \$k+i = Ar* 



y K-l i 3„K i 1„/, _ \„K 



2 r«.-t-i i 2 

The free parameters ^k+Ii^k+z can be easily determined from the first two equations. 
From the final equation we can then compute A. Upon recalling from fl30l) that z can be 
expressed in terms of A, and conversely that A = (z 2 + l)/z + 2, we obtain a polynomial 
equation of degree five for z, 

q(z) := 4z 5 - 12z 4 + 9z 3 - 3z 2 - 4z + 2 = 0. 

Mathematica was unable to factorize q symbolically, hence we computed its roots 
numerically to twenty digits precision. It turns out that q has three real roots and two 
complex roots. The largest real root is at z ~ 2.206272296 which gives the value 

z 2 + 1 

A = + 2 w 4.659525505897. 

z 

The relative errors that we had previously neglected are in fact of order z~ 2K , and 
hence we obtain 

\ K = A, + 0(e~ cK ), where A* w 0.6595 and c w 1.5826. 

This concludes the proof of Theorem |2l 

8.2. Proofs of Lemmas 1141 and 1171 In this appendix, we prove two technical lemmas 
from Section IBTTl Throughout, the iteration matrix G qc \(a) is given by 

G ncl {a):=I-a{A F L)- l Lf, 

where a > and Ap = <fi" F + 4<// 2 ' F > 0. We begin with the proof of Lemma [TH which 
is more straightforward. 



Proof of Lemma 14 Using the basic definition of the operator norm, and the fact that 



Lz = —z'' we obtain 



|Gqci(a)|| ,, 2j00 = max ||(G qcl (a)tt)"|| foo = max II - LG qcl (a)w|| 

II U I Loo =1 I Loo =1 
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We write the operator — LG 



qcl 



a 



-L + -r-IfS as follows: 



[- LG qcl (a)u] i 



u 



if£eC, 



<-f- F (« + <i>WU + K + <+i)) . if ' e A. 

In the continuum region, we simply obtain 



(31) 



[ - LG qc i(a)u] £ = (1 



am a 



for £ G C. 



If t G *4, we manipulate (13T1) . using the definition of Ap = <fi" F + 402 F , which yields 



"f + 2^' F ) 



[ - LG CLd (a)u] £ = 

In summary, we have obtained 

[1 - a]«J', 



S-rh" 

Ap^F 



if £ G c, 
(u'U + u'l +1 ) if£eA. 



— a- 



[ - LG qcl (a)u] 



It is now easy to see that 



< max 



tit 



a .1 



a- 



a 



(1 



)l + 



a 



As a matter of fact, in view of the estimate 



a 



> 11 — a\ — a 



A F 



+ a 



a\ 



the upper bound can be reduced to 

||G ? qc i(«)|| L ( W 2,o 0j W 2,oo) < |l - a(l 



)l + 



a 



(32) 



To show that the bound is attained, we construct a suitable test function. We define 
u via 



sign 



— a 



sign 



(note that G *A for any if > 0) and the remaining values of u" in such a way that 
Yle=-N+i u e = ®- If N > 4, then there exists at least one function u Eli with these 
properties and it attains the bound (I32j) . Thus, the bound in (l32l) is an equality, which 
concludes the proof of the lemma. □ 

Before we prove Lemma [T7] we recall an explicit representation of L _1 L^ cf that was 
useful in our analysis in [TO] . The proof of the following result is completely analogous 
to that of [323 Lemma 14] and is therefore sketched only briefly. It is also convenient for 
the remainder of the section to define the following atomistic and continuum regions 
for the strains: 

A' = { -K + 1,...,K} and C = { - N + 1, . . . , N} \ A'. 
Lemma 21. Let u & U and z = L~ 1 L F ci u, then 



z' e = cr(u')i - cr(n') + <\)" 2F {a_ K {^)h_ K ^ - a K (u')h Ki t), 
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where cr(u'), h±K G M. 2N and cr(V), o>±k{u') G M are defined as follows: 
/ /n = U>e + <t>2F( u e-i + H + u e+i), £G^', 

CTl j ' \(^ + 40' 2 ' F K, <eC, 

! * 

= 2iV = I^V [wjf+i - *4 - U-Jf+1 + u'-k] 7 

e=-N+i 

a_ K (u') = u'_ K+1 - 2u'_ K + u'_ K _ x , a K {u') = u' K+2 - 2u' K+1 + u' K , and 

= [\{ItsK), £=-N+l,...,±K, 
±K ' e H(-l^eK), £ = ±K + 1,...,N. 



Proof. In the notation introduced above, the variational representation of L^ cf from [TO] 
Sec. 3] reads 



{Lf u,v) = (a(u'),v') + 4>2 F [a- K (u')v- 



K - a K (u v K 



Vu,v G U. 



Using the fact that v±n = and *Y^i v 'i = 0' ^ * s eas Y to see that the discrete delta- 
functions appearing in this representation can be rewritten as 

V±K = {h±K,v'). 

Hence, we deduce that the function z = L _1 L^ cf is given by 

(z',v f ) = (L^ u,v) = (er(tt') + (j)2 F [a- K {u)h-K - a K (u')h K },v') \fv G U. 

In particular, it follows that 

z' = a(u') + (j)2 F [a- K {u')h-K - a K (u')h K } + C, 

where C is chosen so that Y^e z i = 0- Since h±K are constructed so that '^2 e h±x,e = 0, 
we only subtract the mean of a(u'). Hence, C = —a(u'), for which the stated formula 
is quickly verified. □ 



Proof of Lemma 11_ Let «eW with < 1. Setting z = G qc \(a)u } and employing 

Lemma EH we obtain 



Uf 



ii f 



A F 



') + 4>2F(®-K(u')h-K,£ ~ a K (u')h Kt i} 



2 ( U K+1 



Ui 



K U -K+l T" "-KJ 

6t- K {u')h- K ^ + a K {u')h K ,£ 



u 



Re + Sp. 



We will estimate the terms Rp and Sp separately. 

To estimate the first term, we distinguish whether £ G C or £ G A'. A quick 
computation shows that Rp = (1 — ot)u' e for £ G C. On the other hand, for £ G A' we 
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have 



Re 



He 



t+i) 



a 



u. 



u 



i+i> 



eA. 



Since H^'H^ < 1, we can thus obtain 



\RA < 



ot\, 



2<j>" 



a 



eA'. 



(33) 



As a matter of fact, these bounds can be attained for certain £, by choosing suitable 
test functions. For example, by choosing u e U with u' N = sign(l — a) we obtain 
Rn = \1 — a\, that is, Rn attains the bound (1331) . By choosing u elA such that 



Ur 



sign 



V-2F 



1 and u\ = sign ( 1 



a 



we obtain that Ri attains the bound ( 1331) . In both cases one needs to choose the 
remaining free u' e so that \u' e \ < 1 and u' e = 0, which guarantees that such functions 
u really exist. This can be done under the conditions imposed on N and K. 

To estimate Si, we note that this term depends only on a small number of strains 
around the interface. We can therefore expand it in terms of these strains and their 
coefficients and then maximize over all possible interface contributions. Thus, we 
rewrite St as follows: 



St 



ce- 



il 



-K-ll 



■h. 



■KM 



U 



u' K [h K ,i - 



-K [2h—K, 

-2h K j 



+ f ] + u'_ 



K+ll 



■h. 



Kj 



L K+l[~^ ,L K/ T jjT U R+2[ 

This expression is maximized by taking u' e to be the sign of the respective coefficient 
(taking into account also the outer coefficient afy^), which yields 



Kj 



+ \h 



|| + 12/1 



Kj 



§1 + 1/1 



KM 



\Si\ < a\^\\\h_ K/ \ + \2h_ K/ + l\ + \h^ + l\ 

= a|^|{|4/i_^ + £| + |4/i^- £ |}. 

The equality of the first and second line holds because the terms ±| do not change the 
signs of the terms inside the bars. Inserting the values for h± K g, we obtain the bound 



15/1 < 



"4 

I A F 



a{A + 2e-4eK)\^\, 



eA, 



and we note that this bound is attained if the values for u' e , I = 
1, K, K + 1, K + 2, are chosen as described above. 

Combining the analyses of the terms Re and St, it follows that 



-K - 1, —K, -K + 



< max 



II - a\ + a4 

1 1 I A F 



|l- a (l-^)|+a(6 + 2e-4£l0|^|} 
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To see that this bound is attained, we note that, under the condition that K > 3 
and N > K + 3, the constructions at the interface to maximize Se and the con- 
structions to maximize Re do not interfere. Moreover, under the additional condition 
N > max(9, K + 3), sufficiently many free strains u' e remain to ensure that Y2i u 'e = 
for a test function u G U, \\u'\\e^> = 1, for which both Rg and Si attain the stated 
bound. That is, we have shown that 

H Gr qd( a 0|| M i,co = max {|! ~ "I + a4 ^4^ 

|l-a(l-^)| + a(6 + 2s-4s^)^} 

=: max{mc(«),m_4(a)}. 

To conclude the proof, we need to evaluate this maximum explicitly. To this end we 
first define a\ = (1 — -^f-) < 1. For < a < «i, we have 

m A {a) = 1 -a + a(4 + 2e-4:£K)\§z\ 
< 1 - a + a4\§f\ =m c {a), 
that is, ||(jqci(Q!)||wi.°° = mc( a )- Conversely, for a > 1, we have 
m A {a) = a(l + (8 + 2s - AeK) 1 -^) - 1 

= m c (a) + a(^ + 2e-4eK)^\ > m c (a), 

that is, || G qc i(«) ||^i,oo = m A (a). Since, in [«i,l], m>c is strictly decreasing and m A is 
strictly increasing, there exists a unique «2 G [or, 1] such that mc(«2) = m^i^) and 
such that the stated formula for ||Gqci(a)||wi.°° holds. A straightforward computation 
yields the value for a.i = a^' 1 ' 00 stated in the lemma. □ 

8.3. Computation of ||G qce || W fc, P . We have computed ||G qce || W fe, P for k = 0,2,p = 
1, 2, oo, from the standard formulas for the operator norm [TTUSH] of the matrix G qcc 
and LGqccL^ 1 with respect to £§. For k = 1 and p = 2, the norm is also easy to obtain 
by solving a generalized eigenvalue problem. 

The cases k = 1 and p — 1, oo are more difficult. In these cases, the operator norm 
of G qce in U l ' p can be estimated in terms of the ^-operator norm of the conjugate 
operator G = / -(L* 8 )-^? : ^ 2N (see LemmaEHfor an analogous definition 

of the conjugate operator L^ nl : ~R 2N — > M. 2N ). It is not difficult to see that ||G qce ||wi.p = 
||G||^ R 2iv for G = I - (Lf^lf : -)• R2JV where we recall that R 2N = g 
j^2N . ^2g<fi£ = 0} (see Lemma [TTl similarly for an analogous definition of the restricted 
conjugate operator L^ nl : M% N — > M% N ), it follows from (|2J) that we have only computed 
||G qce ||^i,p for p = 1, oo up to a factor of 1/2. More precisely, 

HGqceHw 1 '? < IIGjl^P < 2||(j qce ||w 1 >P 

Finally We note that we can obtain L^ cf from the representation given in Lemma [21] 
and that L^ ce can be directly obtained from (jT4{) . 
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